// C port of MOS Bayesian MM
// Transpose b_w to obtain better cache performance
// Make b_c and b_g integers
// Introduce pthreads
// Spawn threads only once for the entire program
// Read data from input file
// Use GMM clusters as initial mu values
// Change update order and initialization for b_c

// include all helper functions
#include <pthread.h>
#include "_bayesian.cc"
#include "_bayesian3.cc"

extern void cudaComputeWPDF(double* answer1, double* answer2, double* answer3, double _newW1, double _newW2, double _newW3, int row, int col
				,double** _bw, int** _bc, int** _bg, double** _wMatrix, int n, int k, int d);
#define NUMTHREADS 2
double L300 = log(1e-300); // log of 1e-300
double XTOL = 1e-4; // tolerance value

// --- Global variables and structures ---
double newW1; // proposed weight
double newW2; // proposed weight
double newW3; // proposed weight
int wrow; // component id of proposed weight
int wcol; // dimension id of proposed weight
double** bwtr; // attribute weights transpose
int** b_c; // active component indicators
int** b_g; // dominant component indicators
double** wMatrix; // against weight matrix
double* __b_alpha; // temporary holder for b_alpha
double** __b_mu; // temporary holder for b_mu
double** __b_mu_inv; // temporary holder for b_mu_inv
double** __b_w; // temporary holder for b_w
int **x; // data points
gsl_rng * r; // the gsl random number generator
struct thread_data {
	int thread_id;
	int n;
	int k;
	int d;
	int n_start;
	int n_end;
	double result1;
	double result2;
	double result3;
}; // struct to be passed to each thread
pthread_t myThreads[NUMTHREADS]; // the threads
struct thread_data myThreadData[NUMTHREADS]; // private data of each thread
int gbl_count; // global variable to count how many slaves have finished their work
int gbl_exit; // flag to tell all workers to exit
enum {FLAG_EXIT, FLAG_W, FLAG_C}; // enumarated types to use for gbl_exit
pthread_mutex_t count_mutex, exit_mutex; // mutex for gbl_count and gbl_exit
pthread_cond_t count_threshold, exit_signal; // conditional variable for signaling count threshold and exit flag

//#include "cudaHost.h"

/*  Helper method to compute the w pdf
 *  newW	proposed weight
 *  row		component id
 *  col 	dimension id 
 *  _bw 	attribute weights
 *  _bc 	active compoment indicators
 *  _bg 	dominant component indicators
 *  _wMatrix	pre-computed matrix containing sum of other weights (n x d )
 *  n		# data points
 *  k 		# clusters
 *  d 		# dimensions
 */
void computeWPDF(double* answer1, double* answer2, double* answer3, double _newW1, double _newW2, double _newW3, int row, int col, double** _bw, 
				 int** _bc, int** _bg, double** _wMatrix, int n, int k, int d, int n_start, int n_end){
	*answer1=0; *answer2=0; *answer3=0;
	double mult1 = 1.0; double mult2 = 1.0; double mult3 = 1.0;
	// create a copy of w(row,:) and update it using proposed weight
	double updateRatio1 = (1 - _newW1)/(1 - _bw[row][col]);
	double updateRatio2 = (1 - _newW2)/(1 - _bw[row][col]);
	double updateRatio3 = (1 - _newW3)/(1 - _bw[row][col]);
	double* _w1 = new double[d]; double* _w2 = new double[d]; double* _w3 = new double[d];
	for (int j=0; j<d; j++){
		_w1[j] = updateRatio1 * _bw[row][j];
		_w2[j] = updateRatio2 * _bw[row][j];
		_w3[j] = updateRatio3 * _bw[row][j];
	}
	_w1[col] = _newW1; _w2[col] = _newW2; _w3[col] = _newW3;
	for (int i2=n_start; i2<=n_end; i2++){ // for each data point
		if (_bc[i2][row]){ // check if component is active
			for (int j=0; j<d; j++){ // for each attribute
				double tempVal1, nextMult1, tempVal2, nextMult2, tempVal3, nextMult3;
				if (_bg[i2][j] == row){ // if the component is dominant
					tempVal1 = (_w1[j] / (_w1[j] + _wMatrix[i2][j])) ;
					tempVal2 = (_w2[j] / (_w2[j] + _wMatrix[i2][j])) ;
					tempVal3 = (_w3[j] / (_w3[j] + _wMatrix[i2][j])) ;
				} else { // if one of the other active components is dominant
					tempVal1 = (_wMatrix[i2][j] / (_w1[j] + _wMatrix[i2][j])) ;
					tempVal2 = (_wMatrix[i2][j] / (_w2[j] + _wMatrix[i2][j])) ;
					tempVal3 = (_wMatrix[i2][j] / (_w3[j] + _wMatrix[i2][j])) ;
				}
				nextMult1 = mult1 * tempVal1; nextMult2 = mult2 * tempVal2; nextMult3 = mult3 * tempVal3;
				if (nextMult1 > 1e-300)
					mult1 = nextMult1;
				else {
					*answer1 += L300;
					mult1 = nextMult1*1e300;
				}
				if (nextMult2 > 1e-300)
					mult2 = nextMult2;
				else {
					*answer2 += L300;
					mult2 = nextMult2*1e300;
				}
				if (nextMult3 > 1e-300)
					mult3 = nextMult3;
				else {
					*answer3 += L300;
					mult3 = nextMult3*1e300;
				}
			}
		}
	}
	*answer1 += log(mult1); *answer2 += log(mult2); *answer3 += log(mult3);
	delete [] _w1; delete [] _w2; delete [] _w3;
}

/*  Helper method to update active component indicators
 *  n		# data points
 *  k 		# clusters
 *  d 		# dimensions
 */
void updateBC(int n, int k, int d, int n_start, int n_end){
	
	for(int i=n_start; i<=n_end; i++){ // for each data point
		
		// mark all dominant components as active, others as inactive
		setValue(0, b_c[i], k-1);
		b_c[i][k-1] = 1; // set default component as active
		for (int j=0; j<d; j++)
			b_c[i][b_g[i][j]] = 1;
		int __bc[k];
		
		// see if remaining compoments could have been active but did not generate any attribute
		//for (int l=0; l<k-1; l++){ // for each passive component
		for (int l=k-2; l>=0; l--){ // for each passive component
			__bc[l] = b_c[i][l]; // save existing value
			//if(b_c[i][l] = 1)
			//continue;
			b_c[i][l] = 0; // mark current component as inactive
			double _bcij_p = log(__b_alpha[l]); // numerator -- component active but did not influence
			double _bcij_q = log(1 - __b_alpha[l]); // denominator -- component not active
			for (int j=0; j<d; j++) { // for each attribute, compute P(xij | cil=0) and P(xij | cil=1)
				double _p_temp = 0; double _q_temp = 0; double _r_temp = 0;
				for (int j2=0; j2<k; j2++){
					_q_temp += b_c[i][j2] * __b_w[j][j2] * (x[i][j] * __b_mu[j][j2] + (1 - x[i][j]) * __b_mu_inv[j][j2]); // P(xij was generated by cij2)
					_r_temp += b_c[i][j2] * __b_w[j][j2] ;
				}
				_p_temp = _q_temp + __b_w[j][l] * (x[i][j]* __b_mu[j][l] + (1 - x[i][j]) * __b_mu_inv[j][l]);
				_bcij_q += log(_q_temp / _r_temp);
				_bcij_p += log(_p_temp / (_r_temp + __b_w[j][l]));
			}
			b_c[i][l] = __bc[l]; // restore old value
			//b_c[i][l] = ( log(gsl_rng_uniform(r)) <= (_bcij_p - logAdd(_bcij_p, _bcij_q)) ); // Pr(c_il = 1) = num / (num + denom)
			__bc[l] = ( log(gsl_rng_uniform(r)) <= (_bcij_p - logAdd(_bcij_p, _bcij_q)) ); // Pr(c_il = 1) = num / (num + denom)
		}
		for (int l=k-2; l>=0; l--)
			b_c[i][l] = __bc[l];
	}
}

/*  Subroutine for the worker thread
 */
void *doWork(void *arg) {
	struct thread_data *myStuff;
	myStuff = (struct thread_data*) arg;
	
	// wait for the signal from master
	// cerr << "Thread #" << myStuff->thread_id << " going to wait state. " << endl;
	pthread_mutex_lock(&exit_mutex);
	pthread_cond_wait(&exit_signal, &exit_mutex);
	while(true){
		// cerr << "Thread #" << myStuff->thread_id << " received wakeup. " << endl;
		if (gbl_exit == FLAG_EXIT){ // check for exit flag
			pthread_mutex_unlock(&exit_mutex);
			// cerr << "Thread #" << myStuff->thread_id << " exiting. " << endl;
			break;
		}
		pthread_mutex_unlock(&exit_mutex);
		
		// start computation
		if (gbl_exit == FLAG_W)
			computeWPDF(&(myStuff->result1), &(myStuff->result2), &(myStuff->result3), newW1, newW2, newW3, wrow, wcol, bwtr, b_c, b_g, wMatrix, myStuff->n, myStuff->k, myStuff->d, myStuff->n_start, myStuff->n_end);
		else if (gbl_exit == FLAG_C)
			updateBC( myStuff->n, myStuff->k, myStuff->d, myStuff->n_start, myStuff->n_end);
		pthread_mutex_lock(&count_mutex); // lock mutex and update count
		gbl_count++;
		if (gbl_count == NUMTHREADS){ // if i am the last slave worker
			pthread_cond_signal(&count_threshold);
			// cerr << "Thread #" << myStuff->thread_id << ": Sending wakeup to master" << endl;
			// cerr << "Thread #" << myStuff->thread_id << " going to wait state. " << endl;
			pthread_mutex_lock(&exit_mutex); // make sure you can get to wait queue
			pthread_mutex_unlock(&count_mutex); // unlock mutex and return to wait state
			pthread_cond_wait(&exit_signal, &exit_mutex);
		} else {
			pthread_mutex_lock(&exit_mutex);
			pthread_mutex_unlock(&count_mutex); // unlock mutex and return to wait state
			// cerr << "Thread #" << myStuff->thread_id << " going to wait state. " << endl;
			// pthread_mutex_lock(&exit_mutex);
			pthread_cond_wait(&exit_signal, &exit_mutex);
		}
	}
	pthread_exit((void*)0);
}

/*  Compute subroutine for the master thread
 */
void computeWithThreads(double* finalSum1, double* finalSum2, double* finalSum3, double _newW1, double _newW2, double _newW3, int _wrow, int _wcol, double _a1, double _b1){
	
	// set up global variables and signal start of computation to workers
	newW1 = _newW1; newW2 = _newW2; newW3 = _newW3; wrow = _wrow; wcol = _wcol;
	pthread_mutex_lock(&count_mutex);
	gbl_count=0; // set worker count
	// pthread_mutex_unlock(&count_mutex);
	pthread_mutex_lock(&exit_mutex);
	gbl_exit = FLAG_W; // set exit flag
	// cerr << "Master: Broadcasting wakeup to all threads." << endl;
	pthread_cond_broadcast(&exit_signal);
	pthread_mutex_unlock(&exit_mutex);
	
	// wait for the signal from workers
	// pthread_mutex_lock(&count_mutex);
	if (gbl_count<NUMTHREADS){ // make sure you will wake up before waiting
		pthread_cond_wait(&count_threshold, &count_mutex);
		// cerr << "Master: Received wakeup." << endl;
	}
	pthread_mutex_unlock(&count_mutex);
	
	//*finalSum1=0.0; *finalSum2=0.0; *finalSum3=0.0;
	*finalSum1 = log(gsl_ran_gamma_pdf(_newW1, _a1, _b1));
	*finalSum2 = log(gsl_ran_gamma_pdf(_newW2, _a1, _b1));
	*finalSum3 = log(gsl_ran_gamma_pdf(_newW3, _a1, _b1));
	for(int tc=0; tc<NUMTHREADS; tc++){
		*finalSum1 += myThreadData[tc].result1;
		*finalSum2 += myThreadData[tc].result2;
		*finalSum3 += myThreadData[tc].result3;
	}
	//	cerr << *finalSum1 << "\t" << *finalSum2 << "\t" << *finalSum3 << endl;
	
}

/*  Update subroutine for the master thread
 */
void updateWithThreads(double* _b_alpha, double** _b_mu, double** _b_mu_inv, double** _b_w, int _n, int _k, int _d){
	
	// set up global variables and signal start of computation to workers
	copy(__b_alpha, _b_alpha, _k);
	copy(__b_mu, _b_mu, _d, _k);
	copy(__b_mu_inv, _b_mu_inv, _d, _k);
	copy(__b_w, _b_w, _d, _k);
	pthread_mutex_lock(&count_mutex);
	gbl_count=0; // set worker count
	// pthread_mutex_unlock(&count_mutex);
	pthread_mutex_lock(&exit_mutex);
	gbl_exit = FLAG_C; // set exit flag
	// cerr << "Master: Broadcasting wakeup to all threads." << endl;
	pthread_cond_broadcast(&exit_signal);
	pthread_mutex_unlock(&exit_mutex);
	
	// wait for the signal from workers
	// pthread_mutex_lock(&count_mutex);
	if (gbl_count<NUMTHREADS){ // make sure you will wake up before waiting
		pthread_cond_wait(&count_threshold, &count_mutex);
		// cerr << "Master: Received wakeup." << endl;
	}
	pthread_mutex_unlock(&count_mutex);
	
}

/*  Estimate beta distribution fit using 3 point values
 */
void estimateBeta(double* _ba, double* _bb, double* _blk, double x1, double x2, double x3, double fw1, double fw2, double fw3){
	double z_a = fw1 - fw2 + log(x1) - log(x2) + log(1-x1) - log(1-x2) ;
	double z_b = fw1 - fw3 + log(x1) - log(x3) + log(1-x1) - log(1-x3) ;
	double z_c = log(x1) - log(x2);
	double z_d = log(x1) - log(x3);
	double z_e = log(1-x1) - log(1-x2);
	double z_f = log(1-x1) - log(1-x3);
	double b_a_num = z_a * z_f - z_b * z_e ;
	double b_b_num = z_b * z_c - z_a * z_d ;
	double denom = z_c * z_f - z_d * z_e ;
	if ((b_a_num > 0 && b_b_num > 0 && denom > 0) || (b_a_num < 0 && b_b_num < 0 && denom < 0)) {
		*_ba = b_a_num / denom; // estimate of a
		*_bb = b_b_num / denom; // estimate of b
		*_blk = (*_ba - 1) * log(x1) + (*_bb - 1) * log(1 - x1) - fw1; // estimate of log(k)
	} else {
		// cout << "Failed Beta: " << b_a_num << "\t" << b_b_num  << "\t" << denom << endl;
		// cout.flush();
		*_ba = log(0); *_bb = log(0); *_blk = log(0);
	}
}

/* Provides three new points to re-estimate beta parameters
 * Points are computed at mean, and one std dev away from mean on both sides
 */
void newBetaPoints(double* _x1, double* _x2, double* _x3, double _ba, double _bb){
	double _bm = _ba / (_ba + _bb); // expected value of the beta
	double _bvar = _ba * _bb / (pow(_ba + _bb, 2) * (_ba + _bb + 1)); // variance of the beta
	*_x2 = _bm;
	*_x3 = _bm + sqrt(_bvar);
	*_x1 = _bm - sqrt(_bvar);
	if (*_x3 > 1)
		*_x3 = 1 - XTOL;
	if (*_x1 < 0)
		*_x1 = XTOL;
}

/*  Main Method
 *  Takes command line arguments for the subsets model
 */
int main(int argc, char* argv[]) {
	
	// --- Start the timer and initialize random number generators ---
	timeval startTime, endTime;
	gettimeofday(&startTime, NULL);
	srand(time(NULL)); // seed srand with start time
	const gsl_rng_type * T; // default init for gsl rng
	gsl_rng_env_setup();
	T = gsl_rng_default;
	r = gsl_rng_alloc(T);
	gsl_rng_set(r,time(NULL)); // seed gsl rng with start time
	cout << setiosflags(ios::fixed) << setprecision(4);
	
	// --- Define which updates are enabled ---
	bool ENABLED_ALPHA = true; // apperance probability updates
	bool ENABLED_C = true; // active components indicators updates
	bool ENABLED_MU = true; // cluster variables updates
	bool ENABLED_W = true; // parameter weights updates
	
	// --- Setup Model Hyper-parameters ---
	double a = 1; // parameter to beta-pdf controlling the alpha parameter
	double b = 1; // parameter to beta-pdf controlling the alpha parameter
	double a1 = 1; // parameter to the gamma-pdf controlling the masks
	double b1 = 1; // parameter to the gamma-pdf controlling the masks
	double a2 = 0.1; // parameter to control the beta-pdf controlling the mean parameters
	double b2 = 0.1; // parameter to control the beta-pdf controlling the mean parameters
	
	// --- Setup Model Parameters ---
	int n=1500; // number of data points
	int k=21; // number of components
	int d=1000; // number of dimensions
	double EPSILON = 1.0/(pow(d,1.0)); // small fixed weight of the default component
	//double EPSILON = XTOL; // small fixed weight of the default component
	cout << "EPSILON = " << EPSILON << endl;
	
	// --- Read input data ---
	allocate(&x, n, d); // allocate memory for data points
	char inFileName[40];
	strcpy(inFileName, "nips1.csv");
	if (readData(n, d, inFileName, x) == -1) {
		cerr << "Call to readData() failed." << endl;
		return -1;
	}
	
	// --- Bayesian Learning Initialization ---
	double *b_alpha; allocate(&b_alpha, k); // learned apperance probabilities
	allocate(&__b_alpha, k);
	//betarnd(b_alpha, k, a, b, r); // init by drawing samples from prior
	setValue(1.0/k, b_alpha, k); 
	b_alpha[k-1]=1; // default component
	
	double **b_mu; allocate(&b_mu, d, k); // learned cluster parameters
	allocate(&__b_mu, d, k); 
	double **b_mu_inv; allocate(&b_mu_inv, d, k); // 1 - b_mu
	allocate(&__b_mu_inv, d, k);
	double* _mu_init; allocate(&_mu_init, d); // holder for init values of mu
	sum(_mu_init, x, n, d, COL); scale(1.0/n, _mu_init, d); // compute mean value of x across each column
	
	for (int i=0; i<k-1; i++){
		int l = gsl_rng_uniform_int(r, n);
		for (int j=0; j<d; j++){
			if (x[l][j] == 1)
				b_mu[j][i] = 0.99;
			else
				b_mu[j][i] = 0.01;
		}
	}
	
	
	for (int i=0; i<d; i++) // for default copy mean(x)
		b_mu[i][k-1] = _mu_init[i];
	//b_mu[i][k-1] = 0.0001;
	copy(b_mu_inv, b_mu, d, k); scale(-1.0, b_mu_inv, d, k); add(1.0, b_mu_inv, d, k); // compute b_mu_inv = 1 - b_mu
	deallocate(&_mu_init);
	
	double **b_w; allocate(&b_w, d, k); // learned parameter weights -- col major order
	allocate(&__b_w, d, k);
	allocate(&bwtr, k, d); // transpose of b_w -- row major order
	setValue(EPSILON, b_w, d, k); // initialize all weights to EPSILON
	setValue(1.0/d, b_w, d, k-1); // initialize weights to 1/d for all components but default
	transpose(bwtr, b_w, d, k);
	
	allocate(&b_c, n, k); // learned active compoments indicators
	
	for (int i=0; i<n; i++) {
		for (int j=0; j<k; j++) {
			b_c[i][j] = (gsl_rng_uniform(r) <= b_alpha[j]) ; 
		}
	}
	
	allocate(&b_g, n, d); // learned dominant component indicators
	for (int i=0; i<n; i++) { // for each data pointer
		for (int j=0; j<d; j++) { // for each attribute
			double f[k];
			if (x[i][j]) {
				for (int l=0; l<k; l++)
					f[l] = b_c[i][l] * b_w[j][l] * b_mu[j][l] ;
			} else {
				for (int l=0; l<k; l++)
					f[l] = b_c[i][l] * b_w[j][l] * b_mu_inv[j][l] ;
			}
			scale(1/sum(f,k),f,k); // weighted multinomial probability for each component
			unsigned int f2[k]; // counts from mnpdf
			gsl_ran_multinomial(r, k, 1, f, f2); // mnpdf with 1 trial
			b_g[i][j] = find(f2,f2+k,1) - f2; // identify the dominant component
		}
	}
	
	// save initialization values
	double *init_alpha; allocate(&init_alpha, k);
	copy(init_alpha, b_alpha, k);
	double **init_mu; allocate(&init_mu, k, d);
	transpose(init_mu, b_mu, d, k);
	double **init_w; allocate(&init_w, k, d);
	transpose(init_w, b_w, d, k);
	int **init_c; allocate(&init_c, n, k);
	copy(init_c, b_c, n, k);
	int **init_g; allocate(&init_g, n, d);
	copy(init_g, b_g, n, d);
	char outFileName[40];
	sprintf(outFileName, "init-c.txt");
	print(init_c, n, k, outFileName);
	sprintf(outFileName, "init-g.txt");
	print(init_g, n, d, outFileName);
	sprintf(outFileName, "init-mu.txt");
	print(init_mu, k, d, outFileName);
	sprintf(outFileName, "init-w.txt");
	print(init_w, k, d, outFileName);
	sprintf(outFileName, "init-alpha.txt");
	print(init_alpha, k, outFileName);
	
	
	// --- Create threads ---
	pthread_mutex_init(&count_mutex, NULL); // init mutex
	pthread_mutex_init(&exit_mutex, NULL);
	pthread_cond_init(&count_threshold, NULL); // init conditional variable
	pthread_cond_init(&exit_signal, NULL);
	double nPerThread = round(((double)n)/NUMTHREADS);
	for (int tc=0; tc<NUMTHREADS; tc++){
		myThreadData[tc].thread_id = tc;
		myThreadData[tc].n = n;
		myThreadData[tc].k = k;
		myThreadData[tc].d = d;
		myThreadData[tc].n_start = (int)nPerThread*tc;
		myThreadData[tc].n_end = (int)(nPerThread*(tc+1) - 1);
		int rc = pthread_create(&myThreads[tc], NULL, doWork, (void*) &myThreadData[tc]);
		if (rc)
			cerr << "pthread_create() failed for thread #" << tc << " with status " << rc << endl;
		else
			cerr << "Created thread #" << tc << endl;
	}
	myThreadData[NUMTHREADS-1].n_end = n-1;
	sleep(2);
	
	// --- Gibbs Sampler ---
	int nitr=3; int mitr=1;
	double **final_alpha; allocate(&final_alpha, nitr, k);
	double ***final_mu; allocate(&final_mu, nitr, k, d);
	double ***final_w; allocate(&final_w, nitr, k, d);
	double *final_ll; allocate(&final_ll, nitr);
	int failed_beta_1 = 0; int failed_beta_2 = 0; int failed_beta_3 = 0;
	for (int gitr=0; gitr<nitr; gitr++) {
		cout << "Iteration #" << gitr << endl;
		
		if (ENABLED_ALPHA){ // update apperance probabilities
			cout << "Updating b_alpha ..." << endl;
			struct timeval st; gettimeofday(&st, NULL);
			int * _bc_sum; allocate(&_bc_sum, k); sum(_bc_sum, b_c, n, k, COL); // sum column-wise b_c
			for (int i=0; i<k-1; i++)
				b_alpha[i] = gsl_ran_beta(r, a + _bc_sum[i], b + n - _bc_sum[i]); // resample from posterior
			deallocate(&_bc_sum);
			struct timeval et; gettimeofday(&et, NULL);
			cout << "Elapsed time = " << timeDiff(et, st) << " seconds" << endl; cout.flush();
		}
		
		if (ENABLED_MU){ // update cluster parameters
			cout << "Updating b_mu ..." << endl;
			struct timeval st; gettimeofday(&st, NULL);
			int ** _bg_sum; allocate(&_bg_sum, d, k); setValue(0, _bg_sum, d, k); // count of dominant component selection
			int ** _x_sum; allocate(&_x_sum, d, k); setValue(0, _x_sum, d, k); // count of generating 1 when dominant
			for (int i=0; i<n; i++){ // for each data point
				for (int j=0; j<d; j++){ // for each attribute
					int tmp_g_ij = b_g[i][j]; // identify dominant component
					_bg_sum[j][tmp_g_ij] += 1 ; // increment corresponding dominant component count
					_x_sum[j][tmp_g_ij] += x[i][j]; // update corresponding x sum
				}
			}
			for (int j=0; j<d; j++){ // for each attribute
				for (int i=0; i<k-1; i++) {// for each component (but default)
					if (_bg_sum[j][i] == 0) // if no data generated then resample from prior
						b_mu[j][i] = gsl_ran_beta(r, a2, b2); // note init_mu is in COL order
					else {// else resample from posterior
						b_mu[j][i] = gsl_ran_beta(r, a2 + _x_sum[j][i], b2 + _bg_sum[j][i] - _x_sum[j][i]);
					}
					b_mu_inv[j][i] = 1 - b_mu[j][i]; // update 1-b_mu
				}
			}
			deallocate(&_bg_sum); deallocate(&_x_sum);
			struct timeval et; gettimeofday(&et, NULL);
			cout << "Elapsed time = " << timeDiff(et, st) << " seconds" << endl; cout.flush();
		}
		
		if (ENABLED_W) { // update parameter weights
			cout << "Updating b_w ..." << endl;
			struct timeval st; gettimeofday(&st, NULL);
			transpose(bwtr, b_w, d, k); // transpose of b_w
			int * _bc_sum = new int[k]; sum(_bc_sum, b_c, n, k, COL); // count of active data points for each component 
			for (int i=0; i<k-1; i++) { // for all components but default
				if (_bc_sum[i] == 0) { // no active data points
					double* _wi = new double[d];
					gammarnd(_wi, d, a1, b1, r); // resample from prior distribution
					scale(1/sum(_wi,d), _wi, d); // normalize
					for (int j=0; j<d; j++)
						bwtr[i][j] = _wi[j]; // update parameter weights
					delete [] _wi;
				} else { // most likely case -- we have active data points
					allocate(&wMatrix, n, d); // against-weight matrix
					for (int l=0; l<n; l++) // setup the matrices
						if (b_c[l][i]){
							b_c[l][i] = 0;
							for (int j=0; j<d; j++){
								wMatrix[l][j] = 0;
								for (int l2=0; l2<k; l2++)
									wMatrix[l][j] += bwtr[l2][j] * b_c[l][l2];
							}
							b_c[l][i] = 1;
						}
					for (int j=0; j<d; j++)// for each attribute
						{
						cout<<"Attribute "<<j<<"Component "<<i<<" ";
						bool FAILED_BETA = false;
						double x1 = 0.9*bwtr[i][j]; // start with a guess for x1,x2,x3
						double x2 = bwtr[i][j];
						double x3 = 1.1*bwtr[i][j];
						if (x3 > 1)
							x3 = 1 - XTOL;
						double fw1, fw2, fw3 = 0;
//arun
//for(int as=0;as<n;as++)
//for(int sa=0;sa<k;sa++)
//b_c[as][i]=1;
						 struct timeval c_start; gettimeofday(&c_start, NULL);
			//			computeWithThreads(&fw1, &fw2, &fw3, x1, x2, x3, i, j, a1, b1); // compute densities at x1, x2, x3
						//Computing using GPU Begins
						cudaComputeWPDF(&fw1, &fw2, &fw3, x1, x2, x3, i, j,b_w, b_c, b_g,wMatrix, n, k, d);
						fw1+=log(gsl_ran_gamma_pdf(x1, a1, b1));
						fw2+=log(gsl_ran_gamma_pdf(x2, a1, b1));
						fw3+=log(gsl_ran_gamma_pdf(x3, a1, b1));
						//Computing using GPU Ends
<<<<<<< .mine
						 struct timeval c_end; gettimeofday(&c_end, NULL);
			//			cout << "Cuda Elapsed time = " << timeDiff(c_end, c_start) << " seconds" << endl; cout.flush();
					//	exit(1);
=======
						 struct timeval c_end; gettimeofday(&c_end, NULL);
			//			cout << "Cuda Elapsed time = " << timeDiff(c_end, c_start) << " seconds" << endl; cout.flush();
						exit(1);
>>>>>>> .r10
						double beta_a1=0; double beta_b1=0; double beta_lk1=0;
						estimateBeta(&beta_a1, &beta_b1, &beta_lk1, x1, x2, x3, fw1, fw2, fw3); // first beta estimate
						if (!isinf(beta_a1)){
							newBetaPoints(&x1, &x2, &x3, beta_a1, beta_b1); // resample x1,x2,x3 using estimated beta
						//	computeWithThreads(&fw1, &fw2, &fw3, x1, x2, x3, i, j, a1, b1); // compute densities at x1, x2, x3
							/*//Computing using GPU Begins
        	                                        cudaComputeWPDF(&fw1, &fw2, &fw3, x1, x2, x3, i, j,b_w, b_c, b_g,wMatrix, n, k, d);
	                                                fw1+=log(gsl_ran_gamma_pdf(x1, a1, b1));
                	                                fw2+=log(gsl_ran_gamma_pdf(x2, a1, b1));
                        	                        fw3+=log(gsl_ran_gamma_pdf(x3, a1, b1));
                                	                // Computing using GPU Ends
*/
							double beta_a2=0; double beta_b2=0; double beta_lk2=0;
							estimateBeta(&beta_a2, &beta_b2, &beta_lk2, x1, x2, x3, fw1, fw2, fw3); // second beta estimate
							if (!isinf(beta_a2)){
								newBetaPoints(&x1, &x2, &x3, beta_a2, beta_b2); // resample x1,x2,x3 using estimated beta
						//		computeWithThreads(&fw1, &fw2, &fw3, x1, x2, x3, i, j, a1, b1); // compute densities at x1, x2, x3
<<<<<<< .mine
							cout<<"There\n";
							//Computing using GPU Begins
                        		                        cudaComputeWPDF(&fw1, &fw2, &fw3, x1, x2, x3, i, j,b_w, b_c, b_g,wMatrix, n, k, d);
                                        		        fw1+=log(gsl_ran_gamma_pdf(x1, a1, b1));
                                                		fw2+=log(gsl_ran_gamma_pdf(x2, a1, b1));
                                             			fw3+=log(gsl_ran_gamma_pdf(x3, a1, b1));
                                               		//Computing using GPU Ends
								cout<<"Here\n";
										double beta_a3=0; double beta_b3=0; double beta_lk3=0;
=======
							/*	//Computing using GPU Begins
                        		                        cudaComputeWPDF(&fw1, &fw2, &fw3, x1, x2, x3, i, j,b_w, b_c, b_g,wMatrix, n, k, d);
                                        		        fw1+=log(gsl_ran_gamma_pdf(x1, a1, b1));
                                                		fw2+=log(gsl_ran_gamma_pdf(x2, a1, b1));
                                             			fw3+=log(gsl_ran_gamma_pdf(x3, a1, b1));
                                               		*/	//Computing using GPU Ends

								double beta_a3=0; double beta_b3=0; double beta_lk3=0;
>>>>>>> .r10
								estimateBeta(&beta_a3, &beta_b3, &beta_lk3, x1, x2, x3, fw1, fw2, fw3); // third beta estimate
								if (!isinf(beta_a3)){
									double _bwij = gsl_ran_beta(r, beta_a3, beta_b3); // sample weight parameter using beta estimate
									double updateRatio = (1 - _bwij)/(1 - bwtr[i][j]); // update all other weights for the component
									for (int j2=0; j2<d; j2++)
										bwtr[i][j2] *= updateRatio;
									bwtr[i][j] = _bwij;
								} else {
									FAILED_BETA = true;
									failed_beta_3++;
								}
							} else {
								FAILED_BETA = true;
								failed_beta_2++;
							}
						} else {
							FAILED_BETA = true;
							failed_beta_1++;
						}
						if (FAILED_BETA){
							// cout << "Failed Beta for w[" << i << "][" << j << "]. Don't do anything for now." << endl;
							// cout.flush();
						}
					}
					deallocate(&wMatrix);
				}
			}
			transpose(b_w, bwtr, k, d); // copy back to b_w
			delete [] _bc_sum ;
			struct timeval et; gettimeofday(&et, NULL);
			cout << "Elapsed time = " << timeDiff(et, st) << " seconds" << endl; cout.flush();
		}
		
		if (ENABLED_C) { // update component indicators
			cout << "Updating b_c ..." << endl;
			struct timeval st; gettimeofday(&st, NULL);
			
			updateWithThreads(b_alpha, b_mu, b_mu_inv, b_w, n, k, d); // update using threads
			
			struct timeval et; gettimeofday(&et, NULL);
			cout << "Elapsed time = " << timeDiff(et, st) << " seconds" << endl; cout.flush();
		}
		
		// update dominant component indicators
		cout << "Updating b_g ..." << endl;
		struct timeval st; gettimeofday(&st, NULL);
		for (int i=0; i<n; i++) { // for each data point
			for (int j=0; j<d; j++) { // for each attribute
				double f[k];
				if (x[i][j]==1) {
					for (int l=0; l<k; l++)
						f[l] = b_c[i][l] * b_w[j][l] * b_mu[j][l]; // use b_mu is x=1
				} else {
					for (int l=0; l<k; l++)
						f[l] = b_c[i][l] * b_w[j][l] * b_mu_inv[j][l]; // use 1-b_mu if x=0
				}
				scale(1/sum(f,k),f,k); // weighted multinomial probability for each component
				unsigned int f2[k]; // counts from mnpdf
				gsl_ran_multinomial(r, k, 1, f, f2); // mnpdf with 1 trial
				b_g[i][j] = find(f2,f2+k,1) - f2; // identify the dominant component
			}
		}
		struct timeval et; gettimeofday(&et, NULL);
		cout << "Elapsed time = " << timeDiff(et, st) << " seconds" << endl; cout.flush();	
		
		// compute log-likelihood of data given the dominant component indicators
		double b_ll=0;
		for (int i=0; i<n; i++)
			for (int j=0; j<d; j++)
				if (x[i][j])
					b_ll += log(b_mu[j][b_g[i][j]]);
				else
					b_ll += log(b_mu_inv[j][b_g[i][j]]);
		cout << "Log-likelihood of data, given the dominant components = " << b_ll << endl; cout.flush();
		
		// save values of b_alpha, b_mu, b_w, and b_ll
		copy(final_alpha[gitr], b_alpha, k);
		transpose(final_mu[gitr], b_mu, d, k);
		transpose(final_w[gitr], b_w, d, k);
		final_ll[gitr] = b_ll;
		
		// print alpha values
		print(b_alpha,k); cout.flush();
	}
	
	// char outFileName[40];
	double* mean_alpha; allocate(&mean_alpha, k); mean(mean_alpha, mitr, nitr-1, final_alpha, nitr, k);
	sprintf(outFileName, "alpha.txt");
	print(mean_alpha, k, outFileName);
	sprintf(outFileName, "ll.txt");
	print(final_ll, nitr, outFileName);
	double** mean_mu; allocate(&mean_mu, k, d); mean(mean_mu, mitr, nitr-1, final_mu, nitr, k, d);
	sprintf(outFileName, "mu.txt");
	print(mean_mu, k, d, outFileName);
	double** mean_w; allocate(&mean_w, k, d); mean(mean_w, mitr, nitr-1, final_w, nitr, k, d);
	sprintf(outFileName, "w.txt");
	print(mean_w, k, d, outFileName);
	deallocate(&mean_alpha); deallocate(&mean_mu); deallocate(&mean_w);
	
	sprintf(outFileName, "final-c.txt");
	print(b_c, n, k, outFileName);
	sprintf(outFileName, "final-g.txt");
	print(b_g, n, d, outFileName);
	sprintf(outFileName, "final-mu.txt");
	print(b_mu, d, k, outFileName);
	sprintf(outFileName, "final-w.txt");
	print(b_w, d, k, outFileName);
	sprintf(outFileName, "final-alpha.txt");
	print(b_alpha, k, outFileName);
	
	
	// --- Output failed beta information ---
	cout << "First Beta approximation failed " << failed_beta_1 << " times out of " << nitr*(k-1)*d << " attempts." << endl;
	cout << "Second Beta approximation failed " << failed_beta_2 << " times out of " << nitr*(k-1)*d - failed_beta_1 << " attempts." << endl;
	cout << "Third Beta approximation failed " << failed_beta_3 << " times out of " << nitr*(k-1)*d  - failed_beta_1 - failed_beta_2 << " attempts." << endl;
	cout.flush();
	
	
	// --- Clean up worker threads ---
	pthread_mutex_lock(&exit_mutex);
	gbl_exit = FLAG_EXIT; // set exit flag
	cerr << "Master: Broadcasting exit to all threads." << endl;
	pthread_cond_broadcast(&exit_signal);
	pthread_mutex_unlock(&exit_mutex);
	sleep(2);
	
	pthread_mutex_destroy(&count_mutex);
	pthread_cond_destroy(&count_threshold);
	pthread_mutex_destroy(&exit_mutex);
	pthread_cond_destroy(&exit_signal);
	cerr << "Master: Exiting." << endl;
	
	// --- Free up memory ---
	deallocate(&x);
	deallocate(&b_alpha);
	deallocate(&__b_alpha);
	deallocate(&b_mu);
	deallocate(&__b_mu);
	deallocate(&b_mu_inv);
	deallocate(&__b_mu_inv);
	deallocate(&b_w);
	deallocate(&__b_w);
	deallocate(&bwtr); 
	deallocate(&b_c);
	deallocate(&b_g);
	deallocate(&init_alpha);
	deallocate(&init_mu);
	deallocate(&init_w);
	deallocate(&init_c);
	deallocate(&init_g);
	deallocate(&final_alpha);
	deallocate(&final_mu,nitr);
	deallocate(&final_w,nitr);
	deallocate(&final_ll);
	
	gettimeofday(&endTime, NULL);
	cout << "Run time = " << timeDiff(endTime, startTime) << " seconds" << endl; cout.flush();
	gsl_rng_free(r);
	return 0;
	
} // end of main()
